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Abstract 



We review and discuss recent advances in the simulation of bulk critical 
phenomena in model fluids. In particular we emphasise the extensions to 
finite-size scaling theory needed to cope with the lack of symmetry between 
coexisting fluid phases. The consequences of this asymmetry for simulation 
measurements of quantities such as the particle density and the heat capacity 
are pointed out and the relationship to experiment is discussed. A general 
simulation strategy based on the finite-size scaling theory is described and 
its utility illustrated via Monte-Carlo studies of the Lennard-Jones fluid and 
a two-dimensional spin fluid model. Recent applications to critical polymer 
blends and solutions are also briefly reviewed. Finally we consider the outlook 
for future simulation work in the field. 
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I. INTRODUCTION 



That the vapour pressure curve of a simple fluid terminates in a critical point has been 
known since the latter half of the last century when Andrews demonstrated the phenomenon 
of critical opalesence in carbon dioxide |]. Shortly afterwards, van der Waals published a 
mean field type theory of liquid- vapour coexistence that predicted the existence of a critical 
point with divergent compressibility, and formed much of the basis of the understanding of 
the critical region for many years thereafter. Arguably, however, the modern era of critical 
phenomena began in the 1940s, when Guggenheim realised that the coexistence curves of 
many simple fluids are not in fact parabolic as predict by van der Waals theory, but nearly 
cubic. At around the same time, Onsager published his famous solution to the 2D ferromag- 
netic Ising model, which among other things allowed it to be demonstrated that the critical 
exponents are strongly non-classical. It was not until the 1960s, however, that activity in 
the field really began to take off. Accurate experimental studies of magnetic systems on the 
one hand, and series expansion studies of model spin systems on the other, providing vivid 
illustration of the phenomena of scaling and universality, already hinted at previously by 
Guggenheim for fluids 0. These findings heralded the beginning of a prodigious effort to 
understand more deeply the nature of critical phenomena, and in particular of finding ways 
of dealing with the profusion of length scales that characterise the critical region and which 
underpin the singularities in physical observables that are its hallmark. Milestones in this 
quest to date, are the renormalisation group || and conformal invariance theories ||, the 
introduction of which have added substantially to the armoury of the critical point theorist. 

Notwithstanding the great strides made in the understanding of critical phenomena over 
recent years, it is the exception rather than the rule that the available theoretical apparatus 
permits exact calculation of the critical properties of a given model system. Many of the 
systems of interest to physicists have steadfastly resisted efforts to find accurate analytical 
treatments of their critical properties. Indispensable, therefore, in tackling these analytically 
recalcitrant models has been computer simulation. Methods such as finite-size scaling or the 
Monte-Carlo renormalisation group permit one, in principle, to extract accurate estimates for 
bulk critical properties from simulations of finite-sized systems. Until recently, however, such 
simulation studies of critical phenomena were principally confined to simple lattice-based 
magnetic systems like the O(N) or q-state Potts models. More complex systems such as off- 
lattice fluids were deemed to be simply too computationally demanding. Gradually, though, 
and with the advent of new simulation methodologies, and ever more powerful computers, 
this situation is changing and it is now becoming possible to perform detailed simulation 
investigations of critical phenomena in a variety of simple and complex fluid systems. These 
developments are opening up for simulation study, a whole new vista of interesting critical 
and phase coexistence behaviour. Phenomena of topical interest encompass not only the 
standard liquid- vapour or liquid-liquid (consolute) critical phenomena, but also multicritical 
and crossover behaviour in diverse classes of systems ranging from simple atomic fluids to 
electrolytes and polymer blends. 

In this article we review recent progress in the simulation of bulk criticality in model fluid 
systems. We begin by describing recent advances in finite-size scaling (FSS) theory for fluids. 
These advances generalise the FSS techniques previously developed in the magnetic context 
in order to take account of the 'field mixing' phenomenon that manifests the lack of symmetry 
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between coexisting fluid phases. This broken symmetry, which is a fundamental issue in 
the critical behaviour of fluids, is shown to have important consequences for simulation 
measurements of quantities such as the particle density, or the specific heat. The relationship 
to experimental studies of field mixing is also discussed. We then describe a simulation 
methodology based on the insights derived from the FSS theory, and illustrate its utility by 
describing, in some detail, recent studies of critical phenomena in the Lennard- Jones fluid 
and tricritical phenomena in a two-dimensional spin fluid model. Attention is also drawn to 
recent FSS studies of critical behaviour in polymer blends and solutions. Finally we assess 
the prospects for further work in the field and highlight a number of outstanding unresolved 
issues. 



II. FINITE-SIZE SCALING THEORY FOR NEAR-CRITICAL FLUIDS 

In the neighbourhood of a critical point, the correlation length £ grows extremely large 
and often exceeds the linear size L of the simulated system. When this occurs, the singu- 
larities and discontinuities that characterise critical phenomena in the thermodynamic limit 
are smeared out and shifted 0. Unless care is exercised, such finite-size effects can lead to 
serious errors in computer simulation estimates of critical point parameters. 

To deal with these problems, finite-size scaling (FSS) techniques have been developed 
||. Use of FSS methods enable one to extract accurate estimates of infinite- volume quan- 
tities from simulations of finite-size. As its basic tenet, the FSS hypothesis holds that for 
sufficiently large £ and L, the coarse-grained properties of a given near-critical system are 
universal M and depend (up to certain non-universal factors) only on specific combinations 
of L and the relevant scaling fields that measure deviations from criticality ||. Systems hav- 
ing short-ranged interactions and a single component order-parameter (such as many fluids 
and magnets) belong to the Ising universality class ie. their scaling functions and 



critical exponents are identical to those of the simple Ising model. However, qualitatively 
different systems, such as those with multi-component order-parameters may exhibit quite 
different universal behaviour. 

For the purpose of setting out the general FSS concepts, however, we shall remain with 
the Ising universality class, for which the critical behaviour is controlled by two relevant 
scaling fields, which we denote r and h. For the Ising model itself, (or equivalently, the 
ordinary lattice gas) there is a special 'particle-hole' symmetry with respect to change of 
sign of the magnetic field H. This symmetry implies that fluctuations in the order-parameter 
5m = m — m c , and the energy density 5u = u — u c are orthogonal (statistically independent), 
ie. (5m5u) = 0. As a result, the phase boundary of the Ising model coincides with the H = 
axis of the phase diagram and thus it is natural to assign r = T — T c and h = H — H c 0. 

Conjugate to the scaling fields r and h are scaling operators £ and A4. For the Ising 
model, one has £ = u (the energy density) and M. = m (the magnetisation). In the near- 
critical region both £ and M. fluctuate strongly and thus their coarse grained (large length 
scale) properties are expected to exhibit universal scaling behaviour. General finite-size 
scaling arguments [PHI^J6[1 predict that the joint distribution pl(A4,£) exhibits scaling 
behaviour of the form 

p L {M,£) ~ A^Atp Mt£ {A^5M, A+5£, A M h, A £ r) (la) 
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where 



A £ = a s L l ' v A M = a M L d ~^ A M A + M = A £ A+ = L d (lb) 

and 

8M = M- <M> C 58 = £- < £ > c (lc) 

where d is the system dimensionality and v and (3 are the standard critical exponents for 
the correlation length and order-parameter respectively 0. The subscripts c in equations |Tc 
signify that the averages are to be taken at criticality. 

The forms of the L dependences of the operator distributions in equation |T5| arise from 
the relationship between the correlation length £, and the system size. For example, in 
the case of the ordering operator, one has SA4 ~ (r c — t) 13 ~ ^ l3 ^ u . Since, however, at 
criticality, the correlation length is bounded by L, one can simply substitute £ — > L, giving 
5A4 ~ L _/3 / iy . A similar arguement pertains to £. To obtain the dependence of the operator 
distributions on the scaling fields, one assumes that the scaling behaviour depends solely on 
the ratio L/£ 0. The known dependences of £ on r and h then yield the given combinations. 

Given appropriate choices for the non-universal scale factors and a £ (equation |TE]), 
the function Pm,£ is expected to be universal. Precisely at criticality, the scaling fields r 
and h vanish by definition, implying that 

PL (M,£) ~ A-UAtP* M>e (A-U5M,Ap£), (2) 

where PXis(x,y) = pM,£(x,y,0,0) is a function describing the universal and statistically 
scale invariant operator fluctuations characteristic of the critical point. It follows that 
Pm ei x -> u) constitutes a hallmark of a universality class, a fact that (as we shall see) can often 
be exploited to obtain accurate estimates of the critical point parameters of fluid systems. 



Although the generality of the FSS expression, equation la, encompasses all members of 
the Ising universality class, irrespective of the symmetry between the coexisting phases, it 
has long been appreciated |15[ that the form of the fluid scaling fields differ from those of the 



Ising model. Specifically, the reduced symmetry of fluids (ie. the fact that (5p5u) ^ 0, with 
p the density) leads to 'mixed' scaling fields comprising linear- combinations of the reduced 
coupling[](inverse temperature) w and applied (reduced chemical potential) field \x: 

r = w c — w + s(fi — p c ) h = p — p c + r{w c — w), (3) 

where the parameters s and r are non-universal (system-specific) quantities controlling the 
degree of field mixing. In particular, r is identifiable as the limiting critical gradient of 
the coexistence curve in the space of p and w. The role of s is somewhat less tangible; 
it controls the degree to which the chemical potential features in the weak scaling field r. 
The directions of the fluid scaling fields are indicated schematically in the phase diagram of 
figure 0. 



1 For convenience we will adopt the reduced coupling w = J/ksT in this definition of the scaling 
fields, rather than the temperature itself. Both definitions are permissible. 
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As a result of the mixed character of the fluid scaling fields r and h, the respective 
conjugate scaling operators £ and M. are found to comprise linear combinations 0,0] of 
the order-parameter (particle density p) and the energy density u. Specifically, one finds 

M = j^gf [p — su], £ = j^gf [u — rp] . (4) 

We shall henceforth term Ai the ordering operator and £ the energylike operator. For the 
correct choice of the field mixing parameters s and r, the operators satisfy the relation 
(5Ai5£) = 0. The forms of the scaling fields and scaling operators for single component 
systems with and without field mixing are summarised in table |I[ 

The operators Ai and £ are the generalised scaling counterparts of the Ising magnetisa- 
tion and energy. Typically, however, in fluid simulation studies, one is more concerned with 
obtaining estimates for the critical density and energy density. Owing to field mixing, these 
quantities are not expected to exhibit the same FSS behaviour as the Ising magnetisation 
or energy. To elucidate their properties it is expedient to reexpress p and u in terms of the 
scaling operators. Appealing to equation || one finds 

u = £-rM, p — M — s£, (5) 

so that the critical density and energy density distributions are 

Pl(u) = Vl{£ - rM) Pl(p)=Pl(M-s£). (6) 

Now the structure of the scaling form |2| shows that the typical size of the fluctuations in 
the energy-like operator will vary with system size like 5£ ~ £~( d ~V^) ^ where 
a is the specific heat exponent and we have employed the hyperscaling relation dv — 2 — a. 
The typical size of the fluctuations in the ordering operator, on the other hand, vary like 
5A4 ~ L _/3 / iy . From this it is is easy to show [[16] that for a given L, the shape of the energy 
and density distributions can be identified with the distribution of the variable 



X e = aJ^SM cos6 + ag 1 8£ sine, (7) 



with 



t ari e u = _^ L -a- a -P)/v ; tan e = f^-a-a-/?)/* (8) 

ra M a M 

where the subscripts u and p signify that the value of 6 corresponds to the energy density 
and particle density distributions respectively. 

The distributions p(Xq) constitute a spectrum of universal functions (parameterised by 
the value of G) describing the particle density and energy distributions of fluids at finite 



L [|T(|[T7]]. This L dependence arises from the different relative strengths of the critical 
fluctuations in M. and £. Since 1 — a > (5, the critical fluctuations in Ai are stronger than 
those in £, causing the distribution Pl(£) to converge to its average value more rapidly with 
increasing L, than pl{-M-)- Consequently, measurements of the distributions of p and u, 
(each of which represent linear combinations of M. and £), yield L-dependent functional 
forms. 

Figure |2| shows the form of p(Xq) for a representative selection of values of G. The 
distributions were constructed from the joint distribution Pl(tti, u) of the critical Ising model 
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obtained in reference fl6| . Since for the Ising model, A4 — > m and £ — > u, the form of p{Xq) 



is obtained simply by taking linear combinations of the form p^((5mcos0 + <5wsinO). For 
= 0° this yields, trivially, the ordering operator distribution Pl{M), while for = 90° 
the form is that of the energylike operator distribution pi{£). Intermediate between these 
values a range of behaviour is obtained, representing the finite L forms of Pl(p) and Pl{u) 
for fluids. Of course, since pi(M.,£) is universal, then so also are all the finite-size forms of 
p L (p) and p L {u). 

In the limit L — > oo, equation [8] implies that both Q u and P approach zero so that 

Pl(u) = PL (-rM) ~ a M - x rL^ v f M {-a M - X TL^ v 8M) (9a) 
Pl(p) = Pl{M) ~ aM-'L^p^iaM^L^dM) (9b) 

It follows that for any finite s and r, the limiting critical point forms of Pl{p) and Pl{u) 
both match the critical ordering operator distribution p* M (x)= J dyp* M £ (x, y) . This result 
reflects the fact that for sufficiently large L, the critical fluctuations in £ are negligible on 
the scale of those in M.. It should be noted, however, that a quite different state of affairs 
obtains for the critical Ising model where, owing to the absence of field mixing (s = r = 0), 
limL^oopjr^ii) =p L {£). 

This alteration to the asymptotic behaviour of Pl{u) turns out to have important ramifi- 
cations for the manner in which the specific heat is measured. Within the canonical ensemble 
of the Ising model, the specific heat is most readily calculated from the variance of the energy 
fluctuations, which at criticality scales with system size like 

C v = L d ((u 2 } - (u) 2 )/k B T 2 oc L a/U . (10) 

For fluids, however, this is not the correct definition to use since the alteration to the limiting 
form of Pl(u) implies that 

L d ((u 2 } - {u) 2 )/k B T 2 oc V' v . (11) 

which scales asymptotically like the Ising susceptibility (fluid compressibility). To recapture 
the Ising behaviour it is instead necessary to consider the fluctuations of the energy-like 
operator 

C v = L d ({£ 2 } - (£} 2 )/k B T 2 oc L a l\ (12) 

Further practical consequences of field mixing are to be found in the finite-size behaviour 
of the average values of the critical density and energy. Specifically, recall that 

(u) c = (£) c - r(M) e (p) e = (M) c - s(£) c . (13) 

Now, on symmetry grounds, the value of {M) c = I Pl{M)<IM. is independent of L. However, 
no such symmetry condition pertains to Pl(£), whose average value (£) c = J oI£pl{£) at 
criticality varies with system size like 

{£) L C - {£)? ~ L-d"*)/", (14) 

implying the same behaviour for the energy and particle densities: 
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{u) L c - (u)f ~ L"(i-«)/" (15a) 
(P)i - (P)T ~ L-^)/". (15b) 

Thus, even with precise knowledge of the location of the critical point, a single measurement 
cannot yield accurate estimates of the infinite volume values of p c and u c . Instead these 
quantities must be estimated by employing the above scaling relations to extrapolate to the 



thermodynamic limit data from a number of different system sizes [[16 . 

It should also be pointed out, that the finite-size scaling expression, equation 0, is strictly 
only valid in the limit of large L. For smaller system sizes, one anticipates that corrections 
to finite-size scaling associated with non-vanishing values of the irrelevant scaling fields 
will become significant |§. These irrelevant fields take the form a\T e + a^r + . . ., where 
9 is the universal correction to scaling exponent, whose value has been estimated to be 
9 0.54 ±3 for the 3D Ising class Incorporating the least irrelevant of these corrections 
into equation §, one obtains f| 

p L (M,S) ~ A+,A+p* Mi£ (A+,5M,AtS£, ai L- / v ) (16) 

As we shall see in section [IV B| , it is necessary to take account of such corrections to scaling, 
if highly accurate estimates of the critical parameters are desired. 

Finally in this section, we note that the above treatment has been developed (primarily 
for reasons of clarity) in terms of pure fluids of the Ising class. More generally though, the 
critical phenomena for a given fluid system of interest may be different. Thus for example 
the critical behaviour of a binary mixture may be Ising like, but one must allow in general 
for critical fluctuations both in the composition and in the total density. In this case the 
two scaling fields r and h will comprise linear combinations of the three physical field i.e the 
temperature, the chemical potential of one species, and the difference in chemical potential 
between the two species ||19|| . Alternatively one may encounter a different universality class 
altogether, which may have more than just two relevant scaling fields. Examples are ternary 
mixtures of fluids or spin fluids (like that considered in section [IV C| ) which, by virtue of 
their coupled order-parameters, exhibit tricritical phenomena. In such cases the tricritical 
behaviour is characterised by three or more scaling fields, each of which comprises (in general) 
linear combinations of the physical fields. The nature and number of these physical fields 
will, of course, depend upon the specific system under consideration. 



III. THE RELATIONSHIP TO EXPERIMENT 



Experimental studies of critical phenomena in fluids have a long and distinguished his- 
tory, which it would be impossible to review here in any great detail. Instead we refer the 
reader to the recent book by Anisimov PO] for a comprehensive treatment. Here we shall 
focus on just one aspect of fluid criticality, namely the field mixing phenomenon. In so doing 
we shall try to emphasise that simulation can often achieve much more than just mimicking 
experiment. Specifically we argue that by exploiting, rather than simply trying to minimise 
finite-size effects, one can often circumvent the difficulties encountered experimentally. 

Although theoretical predictions concerning the existence of field mixing were confirmed 
experimentally quite some time ago P, 2l]j22 |, it has always proved rather difficult to obtain 
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accurate experimental data for quantities such as the field mixing parameters. The chief 
problem is that of the weakness of the field mixing signature in the experimentally measur- 
able observables. By contrast, a simulation is rather less hamstrung because it permits direct 
access to all physical observables of relevance — notably the exact instantaneous values of the 
particle and energy density, the coupling of which manifests the field mixing phenomenon. 
Moreover, the finite-size of the system, which oftens bedevils computer simulation studies of 
critical phenomena, turns out to be a positive boon for field mixing studies. For finite-sized 
systems, field mixing contributes an L-dependent correction to the limiting forms of Pl{p) 
and Pl{u) (cf. section [TT]). The effect of decreasing the system size is thus to magnify the 
size of the field mixing signature in these distributions. Of course the system size shouldn't 
be too small, or else corrections to finite-size scaling will become significant. In most cases, 
though, it turns out to be possible to attain sufficiently large system sizes so that corrections 
to scaling are minimised, while nevertheless retaining clear field mixing effects in Pl{p) and 
Pl{u)- As a result (and as we will demonstrate in the next section), it is in most cases 
relatively easy to measure accurately the field mixing parameters for the system of interest. 
The situation for experiments is more complicated pippf . There one essentially oper- 



ates in the thermodynamic limit and typically has access only to bare observables such as 
the time-averaged density, or response functions such as the isothermal compressibility or 
isochoric specific heat. For the density, the appearance of the energy operator in equation |5| 
leads to the celebrated singularity in the coexistence diameter |§: 

Pd = ~(Pliq + Pgas) = Pc + CLT + br 1 ^ (17) 

where a is a constant representing the 'rectilinear diameter' and b is closely related to the 
field mixing parameter s. Unfortunately, owing to the smallness of the exponent a ~ 0.11 
for fluids of the Ising universality class, the power of the singular term is close to that of the 
leading regular term. Unless extremely small reduced temperatures are attained, the two 
terms are difficult to distinguish experimentally. Further complications derive from gravity- 
induced density gradients in the experimental sample (which are magnified by the large 



near-critical compressibility), as well as hydrodynamic effects and convection currents [21 



Steps must be taken either to minimise these effects during the experiment, or to correct for 
them in the data analysis. Owing to these difficulties, relatively few experiments [^],^2j have 
been able to convincingly identify the diameter singularity at all, let alone obtain accurate 
values for the field mixing parameter s. However, it is doubtful whether a simulation would 
fare much better in attempting to investigate field mixing in this manner, even without the 
problems of spurious sample inhomogeneities. 

The situation with regard to obtaining field mixing properties from experimental mea- 
surements of the response functions, is also far from straightforward. The generalised com- 
pressibility and specific heat are defined with respect to the scaling field derivatives of the 
operators, whose singular behaviour is: 

Xr=(^) r ~r-~< (18a) 

Xh =(f),~ r " a ' ( 18b ) 
which are the scaling equivalents for fluids of the Ising susceptibility and specific heat 
respectively. For pure fluids, however, the response functions measured in practice are the 
isothermal compressibility and isochoric specific heat: 



S 




(19a) 
(19b) 

Asymptotically, it can be shown that these exhibit the same scaling behaviour as the 
generalised response functions of equation [18]. The only residual effect of field mixing is to 
be found in additional non-asymptotic (sub- dominant) contributions. Specifically, one finds 
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k t = A r" 7 + A\T~ a + ■■■ (20a) 
C v = B r- a + B^- 2a + ■■■ (20b) 

In practice, though, it seems difficult to unambiguously separate the asymptotic and non- 
asymptotic terms from one another and from standard corrections to scaling in the experi- 
mental data. 



IV. SIMULATION STUDIES 

In this section we illustrate the practical application of the finite-size scaling concepts 
outlined in section |I[ We do so by reviewing (in some detail) two recent simulation studies 
that have applied FSS analyses to investigate fluid critical behaviour. The first system 
considered in subsection |IV B| , is the well known Lennard- Jones fluid. We demonstrate 
how measurements of the scaling operator distributions can be analysed within the FSS 
framework to yield highly accurate estimates for the limiting critical point parameters of 
this model. We then turn to a more complicated system, namely a two-dimensional spin 
fluid exhibiting a tricritical point. We show how the FSS concepts can be extended to 
incorporate the added complexity of tricritical phenomena, and detail simulation results for 
the model. Finally in section [IV D| we consider recent applications of FSS techniques to 
simulation studies of critical polymer blends and solutions. 

We begin, however, by briefly discussing some issues relating to the choice of simulation 
ensemble. 



A. Remarks on the choice of simulation ensemble 

The benefits that accrue from a FSS analysis of criticality in fluids are to a large extent 
contingent upon the choice of simulation ensemble to be used. Crucial to the viability of 
the study is a choice of ensemble that adequately provides for the strong near-critical order- 
parameter fluctuations. For liquid-vapour transitions, which are principally characterised 
by density fluctuations, Monte-Carlo simulations [^| within the grand canonical ensemble 
(constant /iV^T-ensemble) p4[ have proved a highly effective approach. The strength of 
the grand canonical ensemble (GCE) stems from the freedom of the total particle number 
to fluctuate p3|. Consequently, the near-critical density fluctuations are observable on the 



largest possible length scale, namely that of the system itself. In principle, though, it 
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is also possible to perform a FSS analysis in the canonical iVV^T-ensemble (in which the 
total density is fixed), by studying density fluctuations within sub-blocks of the total system 
However as described in references pT| , p5f , this approach is plagued with practical 



difficulties and seems to be computationally very intensive. The results emerging from such 
studies have, to date, not been of comparable quality to those obtained from GCE studies. 

Another alternative approach for dealing with density fluctuations, is to employ a con- 
stant pressure (NpT) ensemble ||24]| . Here, the total particle number N is held constant, 
but the density can fluctuate by means of volume transitions. Although FSS usually rest 
upon the idea of comparing the correlation length to the (fixed) linear size of the system, 
its use has recently been extended to the iVpT-ensemble ]28|] by expressing the scaling prop- 
erties not in terms of powers of L, but in powers of N. However, while fully feasible as a 
method for studying liquid- vapour criticality, this approach was found (for various technical 



reasons) to be less efficient than use of the GCE for the study of simple atomic fluids |28 



It should nevertheless prove beneficial when dealing with systems such as polymer solutions 
or electrolyte models for which the GCE particle insertion probability can be prohibitively 
small. 

For liquid-liquid transitions, the near-critical region is characterised by strong concen- 
tration fluctuations. In a simulation one typically considers two species of particles (A and 
B) and common practice is to maintain the total particle number, N = + Nb, constant. 
The analogue of the GCE in this case is the semi-grand canonical ensemble (SGCE), in 
which concentration fluctuations are realised by means of Monte-Carlo moves that swop the 
identity of particles A — > B or B — > A, under the control of a parameter representing the 
chemical potential difference between the two species. The SGCE approach has recently 
been successfully employed in conjunction with a FSS analysis of the ordering operator 
distribution Pl{-M), to study critical phenomena in a symmetric mixtures of square- well 



particles |29| and in the Widom-Rowlinson model |30|| . The approach can also be extended 
to asymmetric mixtures (where field mixing effects must be taken into consideration) and 
even to asymmetric polymer blends as described in section [IV D| . 

Most of the recent work on fluid phase coexistence has been performed using the Gibbs 
Ensemble Monte-Carlo (GEMC) method. This technique, details and applications of which 
are reviewed in reference ||31|| , is very useful for mapping out the phase coexistence envelope 
well away from criticality. However, in its most commonly practiced form, the method 
obtains estimates for the critical parameters by extrapolating a power law to coexistence 
curve data obtained well away from the immediate vicinity of the critical point. Due to 
various pitfalls of this approach p2| , |33| , it has thus far not proved able to provide more than 
rough estimates of critical point parameters. It is also not clear how the GEMC method 
can be combined with FSS techniques. For these reasons we will not consider such studies 
further here. 



B. The Liquid- vapour critical point of the Lennard- Jones fluid 

The Lennard- Jones (LJ) fluid, is arguably the simplest model with the credentials (no- 
tably the symmetry) of a realistic atomic fluid. Its interparticle potential takes the form: 

0(r)=46 LJ [(a/r) 12 -(a/r) 6 ] (21) 
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where €lj is the well-depth (which also defines the reduced temperature T* = ksT/eLj), 
and a is a scale parameter. 

Simulations of the LJ system have recently been carried out using a Metropolis algorithm 
p3| within the grand canonical ensemble ||34|| . As is common practice, the Lennard- Jones 
potential was truncated at a cutoff radius r c = 2.5a, and the potential left unshifted. Five 
system sizes were studied corresponding to L = mr c with m — 3,4, 5, 6, 7. The observables 
recorded in the course of the simulations, were the reduced particle density, p* = L~ d Na 3 and 
the dimensionless energy density u* = L" d (4e£j)~ 1( 3?({r})<7 3 , where $({r}) = J2i<j,j 0(l r i — 
rj|). The joint distribution pi{p*, u*) was accumulated in the form of a histogram. A number 
of preliminary short runs were made for the m = 4 system size in order to locate the liquid- 
vapour coexistence curve using a previous simulation estimate of the critical temperature 



35| . This was achieved by tuning the chemical potential p* until the density distribution 
exhibited a double peaked structure indicative of phase coexistence. Histogram reweighting 
3(| was then used in order to explore the coexistence curve in the neighbourhood of this 



simulation temperature. Such a reweighting scheme, use of which is now standard practice 
in simulation studies of critical phenomena, allows histograms of observables obtained at one 
given T and p to be reweighted to yield estimates corresponding to another T' = T + AT 
and p! = p + Ap. 

To obtain a preliminary estimate of the critical point parameters, the universal matching 
condition for the ordering operator distribution Pl{M) was invoked. As observed in sec- 
tion U, fluid-magnet universality implies that the critical fluid ordering operator distribution 
must match the universal fixed point function p% i (x)=J dyp* M £ (x,y) appropriate 
to the Ising universality class. This latter function is independently known from detailed 
Ising model studies [^]. Thus the apparent critical parameters of the fluid can be estimated 
by tuning the temperature, chemical potential and field mixing parameter s (within the 
reweighting scheme) until Pl{M) collapses onto p* M (x). The result of applying this proce- 
dure for the m = 4 data set is displayed in figure ^ , where the data has been expressed in 
terms of the scaling variable x = aJ^L /3 ^ u (Ai — M. c ) and in line with convention, scaled to 
unit norm and variance. The accord shown corresponds to a choice of the apparent critical 
parameters T*(L) = 1.1853, p*(L) = -2.7843. 

Using this estimate of the critical point, more extensive simulations were performed for 
each of the 5 systems sizes, thus facilitating a full FSS analysis. Reweighting was again 
applied to the resulting histograms in order to effect the matching of Pl{M) to p%i{x), thus 
yielding values of the apparent critical parameters. Interestingly, however, the apparent 
critical parameters determined in this manner were found to be L- dependent. The reason 
for this turns out to be significant contributions to the measured histograms from corrections 
to scaling (cf. equation pl|), manifest as an L-dependent discrepancy between the critical 
operator distributions and their limiting fixed-point forms. In the case of the ordering 
operator distribution the symmetry of the Ising problem implies that the correction 

to scaling function is symmetric in M. — {M). In implementing the matching to p* M (x), one 
thus necessarily introduces an additional symmetric contribution to Pl(-M-) associated with 
a finite value of the scaling field r. This latter contribution has, coincidentally, a form that is 
very similar to that of the correction to scaling function, a fact that makes the cancellation 
of contributions possible. It follows, therefore, that the magnitude of the two contributions 
must be approximately equal. 
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Notwithstanding the complications engendered by corrections to scaling, it is nevertheless 
possible to extract accurate estimates of the infinite-volume critical parameters from the 
measured histograms. The key to accomplishing this is the known scaling behaviour of the 
corrections to scaling, which die away with increasing system size like L~ d l v (cf. equation |16|) . 



Since contributions to Pl{M) from finite values of r, grow with system size like tI}I\ it 
follows that implementation of the matching condition leads to a deviation of the apparent 
critical temperature T*(L) from the true critical temperature T* which behaves like 

T*(oc) - T*(L) oc L- {0+l)/v . (22) 

In figure |4] the apparent critical temperature T*(L) is plotted as a function of L~( e+1 ^ u . 
Clearly the data are indeed very well described by a linear dependence, the least squares ex- 
trapolation of which yields the infinite-volume estimate T* = 1.1876(3). The corresponding 
estimate for the critical chemical potential is p* = —2.778(2). 

The size and character of the contribution of corrections to scaling to the operator distri- 
butions can be seen by plotting their forms at the estimated infinite-volume values of T* and 
//*. Considering first the ordering operator distribution, figure [|(a) shows the critical point 
form of Pl(.M) [expressed in terms of the scaling variable x = aj}L^ u {M - M c )}, for the 
two system sizes m = 4 and m = 7. Also shown is the universal fixed point function p* M (x). 
The corrections to scaling, manifest in the discrepancy between the fluid finite-size data and 
the limiting form, are clearly evident in the figure, especially for the m = 4 system size. 
The corresponding situation for the energy operator distribution Pl{£) = / dM.p L (M.,£) 
is shown in figure ^(b), which is plotted together with the limiting fixed point function 
Pe(y) = I Pm e( x ^y)d x > identifiable as the critical energy distribution of the Ising model, 
and independently known from detailed Ising model studies |16|| . The data have all been 
expressed in terms of the scaling variable y = aJ 1 L d, ~ 1 ' v (£ — £ c ). In this case, one observes 
that the corrections to scaling are noticeably larger than for pl(A4), a fact that presumably 
reflects the relative weakness of critical fluctuations in £ compared to those in M.. 

The values of the critical point field mixing parameters, s and r, are implicit in the 
ordering operator distribution. Their values may be estimated (cf. reference ||16|| ) by treating 
each as a fit parameter which is tuned to best optimise the mapping of the critical operator 
distributions onto their limiting fixed point forms, (cf. figure |J). The resulting estimates 
were, however, found to be slightly L-dependent for the smaller system sizes, an effect that 
presumably stems from corrections to scaling. For the two largest system sizes, though, this 
L-dependence is small and the estimates s = —0.11(1), r = —1.02(1) were obtained. 

Addressing now the critical density and energy distributions, figure |6| shows the mea- 
sured forms of Pl{p) and Pl{u) at the assigned values of the critical parameters. To varying 
degrees, these distributions are asymmetric, a fact which as explained in section [IJ stems 
from field mixing effects. Indeed, the forms of the distributions are consistent with those 
shown in figure |2|. The magnitude of the asymmetry also clearly dies away with increasing 
L, as predicted, although the system sizes are still too small to reveal the limiting behaviour 
Pl{p),Pl{u) — > p* M (x). In figure [7|, the values of (p) c and (u) c corresponding to the dis- 
tributions of figure H are plotted as a function of L~^ l ~ a ^ u (cf. equation [T3|). Although 
no allowances have been made for corrections to scaling (the effects of which are certainly 
much smaller than those of field mixing), the data exhibit, within the uncertainties, a rather 
clear linear dependence. Least-squares fits to the data yield the infinite volume estimates 
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p* = 0.3197(4), u* = -0.187(2). 

The finite-size scaling behaviour of Pl{-M.) and Pl{£) a t the critical point, also serve 
to furnish estimates of the exponent ratios (3/v and \jv characterising the two relevant 
scaling fields h and r. Consideration of the scaling form ^| shows that the typical size of the 
critical fluctuations in the energy-like operator will vary with system size like 5£ ~ L-\ d - l / v ) ^ 
while the typical size of the fluctuations in the ordering operator vary like SA4 ~ L~PI V . 
Comparison of the standard deviation of these distributions as a function of system size 
thus affords estimates of the appropriate exponent ratios. In order to minimise systematic 
errors arising from corrections to scaling, this comparison was performed only for the two 
largest system sizes m = 6 and m = 7. From the measured variance of pi (M. ) for these two 
systems, the estimate (3/v = 0.521(5) was obtained. This compares most favourably with 
the three dimensional (3D) Ising estimate |38] of /3/u = 0.518(7). Given though that no 
allowances were made for corrections to scaling, the quality of this accord is perhaps slightly 
fortuitous. 

Implementing an analogous procedure for Pl{£), yields the estimate \jv — 1.67(7), which 
does not agree to within error with the 3D Ising estimate 1/u = 1.5887(4). In this case, 
however, it is likely that the bulk of the discrepancy is traceable to the high sensitivity 
of Pl{£) with respect to the designation of the field mixing parameter r implicit in the 
definition of £ (cf. equation |4j). In the presence of sizable corrections to scaling, it is 
somewhat difficult to gauge very accurately the infinite volume value of r from the mapping 
of Pl{£) onto p}(y)- Studies of significantly larger system sizes than considered here would 
be necessary to alleviate this problem. 



C. Tricriticality in a two-dimensional spin fluid 

As the second example of fluid critical behaviour we consider the phase behaviour of 
a two-dimensional (2D) 'spin fluid', i.e. a 2D fluid of particles, each of which possesses a 
spin degree of freedom p9|. Physically, such a model is representative of a fluid of magnetic 



molecules lightly absorbed onto a smooth substrate |40| . Owing to the coupling of the mag- 
netism and density fluctuations, however, the phase properties of spin fluids are completely 
different from those of non-magnetic simple fluids |j4l],4"2]| such as the LJ system. 

The simplest spin fluid model is a system of hard discs, each of which carries a spin— | 
magnetic moment. Such a system has a configurational energy given by: 

N N 

*}) = -£ JiTifisiSi + £ Ufa) + HY,s t (23) 

i<j i<j i 

with Si = ±1 , Ufrij) is a hard disc potential with diameter a, and H is an applied magnetic 
field. The distant-dependent spin coupling parameter J{rij) is assigned a square well form: 

J(r) = oo r < a 

J(r) = J a <r < 1.5a 

J(r) = r > 1.5a (24) 

(25) 
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The schematic phase diagram of this model in zero field (H = 0) is depicted in figure |8|. 
For high temperatures, there exists a line of Ising critical points (the so-called 'critical line') 
separating a ferromagnetic fluid phase from a paramagnetic fluid phase. The particle density 
varies continuously across this line. As one follows the critical line to lower temperatures, the 
size of the particle density fluctuations grows progressively. At some point, the fluctuations 
in both the particle density and magnetisation are simultaneously divergent and the system 
is tricritical [|43|| . Lowering the temperature still further results in a phase separation between 
a low density paramagnetic gas and a high density ferromagnetic liquid. For subtricritical 
temperatures, the phase transition between these two phases is first order. 

Owing to the interplay between the density and magnetisation fluctuations, the tricritical 
properties of the spin fluid are expected to differ qualitatively from those on the critical line. 
General universality arguments ]7) predict that for a given spatial dimensionality, fluids with 
short-ranged interactions should exhibit the same tricritical properties as lattice-based spin 
systems such as the spin-1 Blume-Capel model [O. However, since spin fluids possess a 



continuous translational symmetry that lattice models do not, this proposal needs to be 
checked. In addition is if of interest to obtain the field mixing parameters of the model and 
compare them with those for a simple critical fluid. 

Before proceeding however, it is necessary to generalise somewhat the FSS equations of 
section [H]> in order to account for the fact that the tricritical point has not two relevant 
scaling field, but three. In general these scaling fields (which we denote g, A and h') are 
expected to comprise linear combinations of the three thermodynamic fields w = J/kbT, 
/i and H . For the spin fluid model considered here, however, the configurational energy is 
invariant with respect to sign reversal of the spin degrees of freedom and magnetic field. 
This special symmetry implies that the tricritical point lies in the symmetry plane H = 0, 
and that the scaling field h' coincides with the magnetic field H, being orthogonal to the 
— w plane containing the other two scaling fields, g and A. Thus one can write 

h' = H-H t (26a) 
A = (/i - /i t ) + r(w t - w) (26b) 
g = w t - w + s(n - fi t ) (26c) 

where the subscript t signifies tricritical values and the parameters s and r are field mixing 
parameters, controlling the directions of the scaling fields in the fi-w plane. The scaling 
fields g and A are depicted schematically in figure ||(b). One sees that g is tangent to the 
coexistence curve at the tricritical point fil |, so that the field mixing parameter r may be 



identified simply as the limiting tricritical gradient of the coexistence curve. The scaling 
field A, on the other hand, is permitted to take a general direction in the fi-w plane and is 
not constrained to coincide with any special direction of the phase diagram ]15| . 



Conjugate to each of the scaling fields are scaling operators, the form of which are given 

by 

M = m (27a) 

V= — - — [p-su] (27b) 
1 — sr 

£ = — - — [u-rp] (27c) 
1 — sr 
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where m is the magnetisation, p is the particle density and u is the energy density. In 



analogy to equation [Ta|, one can make the following finite-size scaling ansatz for the limiting 



(large L) near-tricritical distribution of pL(p,u,m) 
p L (p, u, m) ~ j——p L (a^L d - yi M, a^L d ' y2 V, a^L^E, ai L yi ti, a 2 L y2 X, a 3 L V3 g) (28) 

where pL is a universal scaling function, the ai are non-universal metric factors and the yi 
are the standard tricritical eigenvalue exponents jyj| , which in analogy to equation la can be 



expressed in terms of ratios of the various tricritical exponents characterising the behaviour 
of p,m and £ in the neighbourhood of the tricritical point. 

Precisely at the tricritical point, the tricritical scaling fields vanish identically and the 



last three arguments of equation |28 can be simply dropped, yielding 



p L (p,u,m) ~ ^-pUa^L^M^^L^V^^L^S), (29) 

where p* L is a universal and scale invariant function characterising the tricritical fixed point. 
In what follows we describe how the tricritical point of the 2D spin fluid model was located 
by means of simulation, and how the proposed universality of equation was tested by 
obtaining the form of p\ and comparing it with that for the tricritical 2D Blume-Capel 
model whose Hamiltonian is given by 

n = - j^^sj + a^sI Si = -1,0,1 (30) 

(u) i 

where A is the so-called crystal field (analogous to the chemical potential). The phase 
diagram of the Blume-Capel model has the same topology as that of figure |8| 

Monte-Carlo simulations of the spin fluid model were performed using a Metropolis 
algorithm within the grand canonical ensemble ||39|| . Three system sizes corresponding to 
L = 18(j, 24<t and 30<r were studied, containing (at liquid- vapour coexistence), average 
particle numbers of (N) « 120, 210 and 330 respectively. During the course of the simulation 
runs, the joint probability distribution pi{p, u, m) was obtained in the form of a histogram. 
The histogram extrapolation technique was used to explore the coexistence region. 

In contrast to the situation for the LJ fluid, no independent knowledge of the tricritical 
scaling functions was available to aid location of the tricritical parameters. It was thus 



necessary to employ the cumulant intersection method to determine T* and p* [13[]. The 
fourth order cumulant ratio Ul is a quantity that characterises the form of a distribution 
and is defined in terms of the fourth and second moments of a given distribution 

U L = 1- -1^. (31) 

3 < m 2 > 2 

The tricritical scale invariance of the distributions pl(T>),Pl(M) and Pl(£), (as expressed 
by equation ^), implies that at the tricritical point (and modulo corrections to scaling), the 
cumulant values for all system sizes should be equal. The tricritical parameters can thus be 
found by measuring Ul for a number of temperatures and system sizes along the first order 
line, according to the prescription outlined below. Precisely at the tricritical temperature, 
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the curves of Ul corresponding to the various system sizes are expected to intersect one 
another at a single common point. 

Initially, a very approximate estimate of the location of the tricritical point was inferred 
from a number of short runs in which a temperature was chosen and the chemical potential 
tuned while observing the density. Long runs were then carried out using this estimate for 
each of the three system sizes. From these runs the first order line and its analytic extension 
in the p-T plane were identified. The measured values of lf£ along this coexistence line 
are shown in figure |9] for the each system size. Clearly the curves of figure ^| have a single 
well-defined intersection point, from which the tricritical temperature may be estimated as 
being fc B T t /J = 0.581(1). The associated estimate for the chemical potential is p t /k B T = 
—1.916(2). The tricritical particle and energy densities were estimated to be p t = 0.374(1) 
and Ut = —0.778 respectively. The average magnetisation is of course strictly zero on 
symmetry grounds. Typical near-tricritical configurations for the L = 30a system are shown 
in figure [T^. They clearly show the coupling of the density and magnetisation fluctuations, 
with a disordered spins spin configuration at low density, and an ordered spin configuration 
at high density. 

Figure [11] shows the forms of the operator distributions Pl(M), PiXP) and Pl(£) corre- 
sponding to the designated values of the tricritical parameters. Also included in figure O are 
the scaled tricritical operator distributions obtained in a separate study of the 2D Blume- 
Capel model |3j|. Clearly in each instance and for each system size, the scaled operator 
distributions collapse extremely well onto one another as well as onto those of the tricritical 
Blume-Capel model. This data collapse is perhaps the most stringent test of universality, 
and transcends critical exponent values, although these too were obtained in the study of 



reference |39| and agree with those of the 2D Blume-Capel model. There can thus be little 
doubt that despite their very different microscopic character, the spin fluid and Blume-Capel 
models do indeed share a common fixed point. 

The values of the field mixing parameters r and s, are implicit in the forms of Pl(£) and 



Pl(T3) shown in figure |IT], from which one finds |39| r = —2.82, s = —0.013. Intriguingly, 
this value of s is at least an order of magnitude smaller than that measured at the critical 
point of the 2D Lennard- Jones fluid [|T2]] . This smallness implies that the scaling field A 
almost coincides with the p axis of the phase diagram. 

With regard to the forms of the tricritical operator distributions, it is found that p\{£ ) is 
(to within the precision of the measurements) essentially Gaussian, implying that the tricrit- 
ical fluctuations in S are extremely weak. The tricritical form of p\{M) on the other hand is 
three-peaked, in stark contrast to the situation in the 2D Ising model, where the critical mag- 
netisation distribution is strongly double-peaked ]T3",i5l|. This three-peaked structure reflects 



the additional coupling between the magnetisation and the density fluctuations. Specifically, 
the central peak corresponds to fluctuation to small density, which are accompanied by an 
overall reduction in the magnitude of the magnetisation (cf. figure [10]). Were one, however, 
to depart from the tricritical point along the critical line, these density fluctuations would 
gradually die out and a crossover to a magnetisation distribution having the double-peaked 
Ising form would occur. 

Finally in this subsection, we mention that recent theoretical and simulation studies 
have also been carried out on a 3D spin fluid model having Heisenberg instead of Ising 
spins p4P>||ET| . GEMC simulations were carried out to determine the liquid- vapour phase 



16 



envelope, while the behaviour at points on the critical line of magnetic transitions was 
investigated using a FSS analysis in the canonical ensemble. Interestingly the results were 
not in full accord with the known properties of the lattice Heisenberg model, indicating a 
possible failure of universality. Further work is clearly called for in order to clarify this issue. 



D. Criticality in polymeric fluids 

Simulations of polymer systems are considerable more exacting in computational terms 
than those of simple liquid or magnetic systems. The difficulties stem from the problems 
of dealing with the extended physical structure of polymers, which gives rise to extremely 
slow diffusion rates, manifest in protracted correlation times. In order to ameliorate these 
difficulties, coarse-grained lattice-based polymer models such as the bond fluctuation model 
(BFM) [Q are often employed. Within the framework of the BFM each monomer occupies 
a whole unit cell of a 3D periodic simple cubic lattice and neighbouring monomers along the 
polymer chains are connected via one of 108 possible bond vectors, providing for 5 different 
bond lengths and 87 different bond angles. The model sacrifices chemically realistic detail 
for greater computational tractability, while nevertheless retaining the essential qualitative 
features of polymeric systems, namely chain connectivity and excluded volume. For studies 
of critical phenomena, this neglect of chemical detail is not expected to bear on the universal 
scaling properties. 

Most studies of criticality in polymeric systems have focussed on symmetric blends i.e. 
two species having identical chains length Na = Nb- SGCE simulations (cf section [IV A|) 



of simple lattice walks were studied by Sariban and Binder H3] who accurately located the 
critical temperature by using the cumulant intersection method (cf. Section [IV Q) . A similar 
approach was applied by Deutsch and Binder PJlHJ to the BFM, this time making extensive 



use of histogram reweighting to map the phase diagram as a function of chain length. They 
confirmed the Ising character of the critical point and reported tentative evidence for a 
crossover from Ising to mean field behaviour away from the critical point |52] . 



Work on strongly asymmetric polymer blends Na ^ Nb, is technically more complicated 
than for symmetric mixtures, and has only become feasible very recently. The difficulties 
stem from excluded volume effects, which prevent one simply removing a chain of one species 
and replacing it with another of unequal size. For fairly short chains (N < 30) this problem 
can be tackled using a novel simulation algorithm due to Miiller and Binder |53|, which 



operates by making MC moves that break an A-chain into k B-chains or joins k B-chains to 
form a single A-chain. This algorithm was recently also employed by Miiller and Wilding 



53 1 to study the critical behaviour of a polymer blend within the BFM, having Na = 
10, Nb = 30. The observables sampled were the concentration of pa of A— type chains (the 
concentration of B chains being fixed by stipulating NaPa + NbPb = constant) and the 



energy density. FSS methods similar to those described in sections [IV B| and |1V (J| were 



employed to accurately locate the critical point and coexistence curve as a function of the 
temperature and chemical potential difference between the two species. The Ising character 
of the critical point was clearly demonstrated and the field mixing parameters obtained. The 
scaling of the critical temperature with chain length asymmetry was further studied using 



FSS by Miiller and Binder [55], who found: 



17 



T c oc N A N B /(N)l 2 + < /2 ) 2 , (32) 

in agreement with Flory-Huggins mean-field theory. For short chains, however, the mean 
field theory was found to overestimate the critical temperature by about 25% due to the 
failure to treat properly the critical fluctuations. 



Very recently, Miiller and Schick |56| have studied a ternary polymer blend comprising 
two incompatible homopolymers and a symmetric diblock copolymer. Such a system has a 
phase diagram somewhat similar to that of the spin fluid considered in section [IV C| , with a 
line of second order transitions (corresponding to criticality between the two homopolymer 
species) ending in a tricritical point below which there is three-phase coexistence between 
two homopolymer rich phases and a spatial structured copolymer-rich phase. The SGCE 
simulation approach was again employed and the Ising character of the critical line confirmed 



using FSS methods similar to those described in the section |V_B[ The approximate location 
of the tricritical point was also determined, although the accessible range of chain lengths 
and system sizes were not sufficiently large to allow an unambiguous characterisation of its 
nature. 

Finally, studies of the polymer- solvent critical point have recently been reported by Wild- 
ing et al ||571| . A biased chain insertion scheme was employed ||58|| , allowing a GCE algorithm 



to be implemented for the BFM . Chains up to length N = 60 monomers were studied and a 
basic FSS study carried out in order to determine the chain length dependence of the criti- 
cal temperature and volume fraction. For each chain length investigated, the critical point 
parameters were determined by matching the ordering operator distribution function to its 
universal fixed-point Ising form. Histogram reweighting methods were employed to increase 
the efficiency of this procedure. The results indicate that the scaling of the critical temper- 
ature with chain length is relatively well described by Flory theory, i.e. — T c ~ iV~ a5 , 
where is the temperature at which the chains behave ideally. The critical volume fraction, 
on the other hand, was found to scale like <p c ~ iV~ a37 , in clear disagreement with the Flory 
theory prediction <p c ~ iV~ a5 , but in good agreement with experiment [1511. Measurements 



of the chain length dependence of the end-to-end distance suggested that the chains are not 
collapsed at the critical point. 



V. DISCUSSION AND OUTLOOK 

In this review we have attempted to show that the finite-size scaling techniques previ- 
ously deployed with such success for lattice spin models, can also be extended to permit 
accurate study of critical behaviour in continuum fluids and polymeric systems. In addition 
to illustrating how one can tackle the critical region in practice, the applications reviewed 
were intended to convey a flavour of the types of systems and phenomena that have been 
studied. Needless to say, however, there are a whole host of other fluid systems which have 
not yet been studied in detail by simulation, and whose critical phenomena raise a number 
of very interesting questions. In this section we attempt to highlight a small selection of 
them. 



One issue that is currently generating great theoretical and experimental interest [60-64 



is the thorny question of the nature of criticality in ionic fluids. Although the Ising character 
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of simple fluids with short range interactions is well established, the theoretical understand- 
ing of ionic criticality is much less developed. Owing to their long ranged coulomb interac- 
tions, ionic fluids are not expected a-priori to belong to the Ising universality class. However 
the possibility of an effective short-ranged interaction engendered by charge-screening, does 
not preclude the possibility of asymptotic Ising behaviour. Indeed, experimentally, examples 
have been found of ionic fluids exhibiting pure classical (mean-field) critical behaviour, pure 
Ising behaviour, and in some cases a crossover from mean field behaviour to Ising behaviour 
as the critical point is approached p0| . It is thus of considerable interest to determine the 



asymptotic universal behaviour, as well as the physical factors controlling the size of the 
reduced crossover temperature. Simulation work on this problem to date has concentrated 
on the prototype model for an ionic fluid, namely the restricted primitive model (RPM) 
electrolyte, which comprises a fluid of hard spheres, half of which carry a charge +q and 
half of which carry a charge —q. Unfortunately, simulations of this model are hampered 
by slow equilibration and the need to deal with the long range coulombic interactions p5 



Nevertheless very recently, a FSS simulation study of the RPM has been reported, which 
applied the techniques of section [IV BL and which appears to show that the critical point 



is indeed Ising like in character ||66|| . It remains to be seen though, to what extent this 
conclusion extends to other less artificial electrolyte models, and whether simulation can 
shed light on the system-specific factors controlling the disparaties in crossover behaviour 
observed in real systems. 

Somewhat similar questions to those concerning ionic fluids, have also recently been 
raised in the context of simple fluids with variable interparticle interaction range. Gibbs 
ensemble Monte Carlo studies of the coexistence curve properties of the square-well fluid 
|67j seem to suggest that while for short-ranged potentials the shape of the coexistence 



curve is well described by an Ising critical exponent j3 ~ 0.324, the coexistence curve for 
longer ranged potentials is near-parabolic in shape implying a classical (mean-field) exponent 
P ~ 0.5. This finding could doubtless be investigated in greater detail using FSS techniques. 
If a crossover does indeed occur, it would be useful to try to formulate a Ginzburg criterion to 
describe its location as a function of the interparticle interaction parameters. Experimental 
and theoretical results concerning crossover from Ising to mean field behaviour in insulating 



fluids and fluid mixtures have also recently been discussed in the literature [68,20 



Another matter of long-standing interest is the question of the physical factors governing 
the size of the field mixing effect in fluids, which seem to be strongly system specific. Thus, 
for instance, the measured size of the coexistence curve diameter singularity (cf. subsec- 
tion p|) in molten metals is typically much greater than that seen in insulating fluids. This 
has prompted the suggestion |)0 that the size of the field mixing effect is mediated by the 
magnitude of three-body interactions. Clearly there is scope for simulation to attempt to 
corroborate this proposal, possibly by introducing many body interactions into a simple fluid 
model and observing the values of the field mixing parameters as their strength is tuned. 

On a slightly different note to those emphasised in this review, we remark that very little 
has been done to study interfacial effect close to the critical point of realistic fluid models. 
Thus there are, to date, no accurate simulation studies of the surface tension critical exponent 
or associated amplitudes ]70]. One might conceivably investigate this matter by studying the 



temperature dependence of the ordering operator distribution Pl{M) at coexistence, using 
preweighting techniques to overcome the free energy barrier, as has recently been done for 
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the Ising model and the Lennard- Jones fluid [fn] ,[72||3-i 



Turning finally to polymer systems, we note that recent advances in simulation tech- 
niques, such as biased growth techniques and chain breaking algorithms contribute signif- 
icantly to our ability to deal with critical phenomena in polymer systems. Nevertheless, 
the studies reviewed in subsection [IV D| were still limited to relatively small chain lengths 
and system sizes and it is questionable whether the asymptotic (large chain length) scaling 
properties are really being probed. While growing computational power will help somewhat 
to alleviate this problem, further algorithmic improvements are clearly called for in order to 
deal with even longer chains. 
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TABLES 





Ising model 


Pure fluids 
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T-T c 


w c - w + s(p - He) 
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H-H c 


fi- fx c + r(w c - w) 
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± } sr [u rp] 


M 


m 


i_a r [P su \ 



TABLE I. The forms of the relevant scaling fields and scaling operators for the Ising magnet 
and for pure fluids. 
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FIG. 1. Schematic representation of the liquid- vapour coexistence curve showing the directions 
of the relevant scaling scaling fields. The angles ip\ and tp2 are related to the field-mixing parameters 
s and r (equation ||) by r = — tan ipi and s = tan if)2 ■ 
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FIG. 2. Selections from the universal finite-size spectrum of critical density and energy density 
distributions of fluids. The distributions were obtained according to the procedure described in 
the text. Following convention, the values of the non-universal scale factors 1 and aj^ have been 
chosen to ensure that the distributions have unit variance. From ref. 1TF 
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FIG. 3. The measured form of the ordering operator distribution pl(-M) for the m = 4 system 
size at the apparent critical parameters T* = 1.1853,//* = —2.7843. Also shown for comparison 
is the universal fixed point ordering operator distribution p* M (x). The data has been expressed 
in terms of the scaling variable x = aJ^L@/ u (A4 — Ai c ), with the value of the non-universal scale 
factor aj^ chosen so that the distributions have unit variance. Statistical errors do not exceed the 



symbol sizes. From ref. [34] 
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FIG. 4. The apparent reduced critical temperature, (as denned by the matching condition 
described in the text), plotted as a function of L - ^ 1 )/", with 9 = 0.54 and u = 0.629 ||,||]. The 
extrapolation of the least squares fit to infinite volume yields the estimate T* = 1.1876(3). From 
ref. @ 



28 



FIG. 5. (a) The ordering operator distribution pl{M) for the two system sizes m = 4 and 
m = 7 at the assigned critical parameters T*,fj,*, expressed as function of the scaling variable 
x = aJ^L@/ u (A4 — A4 C ). Also shown (solid line) is the universal fixed point ordering operator 
distribution p* M {x). (b) The energy operator distribution Pl{£) for the two system sizes m = 4 
and m = 7 at T* , /j,* , expressed as a function of the scaling variable y = a^ 1 L d ~ l / u {E — £ c ). Also 
shown (solid line) is the universal fixed point energy operator distribution p^(y) ■ In both cases 
the values of the non-universal scale factors a^ 1 or aj^ have been chosen to yield unit variance. 
From ref. [p4| 
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FIG. 6. (a) The density distribution at T* , fi* for the system sizes m = 4 — 7, (b) The corre- 
sponding energy density distributions. The lines are merely guides to the eye. Statistical errors do 
not exceed the symbol sizes. From ref. [34] 
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FIG. 7. (a) The measured average density (p) c (L) at the designated critical point, expressed 
as a function of L~( 1 ~ a ^ u . The least-squares fit yields an infinite volume estimate p c = 0.3197(4). 
(b) The measured average energy density {u) c {L) at the critical point, expressed as a func- 
tion of L~( 1 ~ a ^ u . Extrapolation of the least-squares fit to infinite volume yields the estimate 
u c = -0.187(2). In both cases \jv = 1.5887 was taken @. From ref. @ 
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FIG. 8. (a) Schematic phase diagram of the spin fluid in the T—p plane, (b) Schematic phase 
diagram in the fj,—T plane showing the directions of the relevant scaling fields g and A. From ref. 



35 



0.61 



0.59 
0.58 



1 1 


1 1 1 1 

L=18o 




L=24o 




L=30c 


1 1 


ii 



0.56 
0.55 
0.54 
0.53 

0.574 0.576 0.578 0.580 0.582 0.584 0.586 0.588 

K B T/J 

FIG. 9. The measured cumulant ratios Ul for the 2D spin fluid model along the first order line 
and its analytic extension determined according to the procedure described in the text. From ref. 
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FIG. 10. Typical particle/spin configurations of the L ■ 
values of +1 are denoted by filled circles, and spin values 



30<7 spin fluid near tricriticality. Spins 
-1 by unfilled circles. From ref. 
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FIG. 11. The scaling operator distributions for the 2D spin fluid at the designated tricritical 
parameters for each of the three system sizes L = 18a, 24a, 30a. (a) p^(A4), (b), p^iT 1 ) (c) 
p* L {6). Also shown for comparison are the corresponding distribution measured for the tricritical 
L = 40 2D Blume-Capel model. All distributions are expressed in terms of the scaling variable 
a^ 1 L d ~ Vi (O — O c ) and are scaled to unit norm and variance. Statistical errors do not exceed the 



symbol sizes. From ref. [3£] 



39 



